This is the master document organizing all MFC2 data analyses! Think of this as a data summary + table of contents. Specific analyses will be conducted and figured generated within linked RMarkdown files, acting like chapters. With this file, you can easily locate specific analyses of interest.
MFC2 = “Maternal-Fetal Conflict 2” = 2nd pilot study investigating the role of the gut microbiome in mediating maternal-fetal conflict over gestational energy allocation.
We investigated two main questions in this pilot:
How do placental lactogen (PL) and placental growth hormone (PGH) alter the composition of the gut microbiome?
To what extent are the metabolic effects of each hormone dependent on the presence of a gut microbiome?
Here’s an unpolished schematic overview of the study design:
Here’s a schematic overview of the various data types collected and
analyses that were run:
Note: This was probably the least interesting data collected.
Our endpoint phenotypic data consists of the following measurements:
Here’s an interactive table to take a look at the raw endpoint data:
Here’s a summary table of the three endpoint phenotypic parameters across the six groups:
For phenotypic data, our analyses will be broken down as follows:
Before examining each individual metric, I was curious what the
relationship is between colon and SI lengths. Do mice with larger SIs
also tend to have larger colons?
[1]
“The linear model is significant. p = 0.00756 and R^2 = 0.168 .”
For
the ABX mice, The linear model is not significant. p = 0.105 and R^2 =
0.103 . For the Non-ABX mice, The linear model is not significant. p =
0.92 and R^2 = -0.0618 .
Those results were not particularly interesting.
First let’s just take a look at the data across all groups:
NOTES:
Outliers detected and removed: 2 outlier(s).
Takeaways: With outliers excluded, there is a signif difference in SI length between PL and saline-exposed mice with a conventional gut microbiome, but not between those who received antibiotics. No other within-abx-group comparisons reached significance. When the data is left as-is, there are no significant differences between hormone groups.
Now let’s rearrange the data to be faceted by hormone group,
comparing abx and non abx mice.
Takeaways: As expected, mice exposed to antibiotics have longer
small intestines than those with a complete conventional gut microbiome.
However, this difference only reaches significance in Saline and PGH
mice. In PL mice, we approach but don’t quite hit the 0.05 threshold.
This is true regardless of the inclusion or exclusion of outliers.
Takeaways: No observed difference in colon length across
hormone groups. I didn’t bother to remove outliers from this data
because it seemed clear to me that it wouldn’t make any difference.
Now let’s rearrange to compare each hormone group by abx status:
Takeaways: Colon lengths were generally longer in mice exposed to antibiotics in all hormone exposure groups, but this only reached significance in the PL-exposed mice. Note: I reran this excluding the one saline non-abx outlier, but it made no difference. The range of colon lengths within each group was also quite large, especially for the abx-exposed mice.
Takeaways: No observed difference in cecal mass across hormone
groups. I didn’t bother to remove outliers from this data because it
seemed clear to me that it wouldn’t make any difference.
Takeaways: As expected, regardless of hormone exposure, mice
exposed to antibiotics have much larger cecums than those with a
conventional gut microbiome.
Here’s a simple review of the difference between what OGTTs and ITTs tell us:
OGTT: Measures how efficiently an organism processes glucose, assessing glucose tolerance and the overall function of glucose regulatory mechanisms. A mouse’s response to an OGTT depends on both the ability of the pancreas to secrete insulin in response to glucose and the effectiveness of tissues (e.g., muscle, liver, adipose tissue) in taking up glucose out of circulation. Poor glucose processing during an OGTT suggests glucose intolerance.
ITT: Measures how effectively insulin lowers blood glucose levels (i.e. the sensitivity / responsiveness of tissues to insulin). Reduced glucose clearance during an ITT suggests insulin resistance. Together, OGTTs and ITTs provide us with a complementary look at glucose homeostasis.
AUC is a quantitative measure of the total blood glucose exposure over time during an OGTT or ITT. It is calculated by summing the area under the glucose curve generated by plotting glucose levels against time (i.e. integration). AUC allows for easy comparison between groups or treatments by summarizing complex time-series data into a single, interpretable metric.
NOTES:
Two glucose measurements were taken for each mouse at each time point. If the 2 measurements were >25 apart, a 3rd measurement was taken. Mean blood glucose (of the 2-3 measurements) was calculated. Measurements >1 SD away from the mean were removed, and the mean was re-calculated.
While OGTTs have a nearly 100% successful execution rate, ITTs occasionally are unsuccessful (likely due to human error in injection). For this experiment, I set a cutoff for a successful ITT as a >15% drop in blood glucose levels from fasting (time 0) to time 15. Data from tests that failed to meet this threshold were not included in the final data files. Unfortunately, this means that some mice do not have both baseline and endpoint ITT measurements.
While less relevant to the hypotheses of this project, it’s notable that there are differences between responses to the tolerance tests in abx and non-abx mice prior to any hormone injections (at baseline).
Here are two visualizations of those baseline differences. The first
shows the mean response with standard error bars, and the second shows
each individual mouse’s response with the group means overlaid:
Below is the same two visualization styles, but for ITT responses:
Now let’s look at the AUCs for the baseline OGTT and ITTs in ABX and Non-ABX mice:
Is there a significant difference at baseline between ABX and non-ABX OGTT AUCs? [1] “The t-test result is significant. p = 3.17e-06 .”
Is there a significant difference at baseline between ABX and non-ABX ITT AUCs?
[1] “The t-test result is not significant. p = 0.962 .”
Key takeaway: We are observing discordance between the OGTT and ITT data. The OGTT data shows greater glucose tolerance (indicating lower insulin resistance) in the ABX mice (as expected), while in the ITT data, the Non-ABX mice exhibit greater glucose drops (potentially indicating lower insulin resistance). This apparent contradiction may be resolved by looking at both in terms of % change in addition to absolute change:
NOTE: I’m unsure if it would make sense to calculate AUCs for the %-based graphs. There’s certainly an argument for it, but I don’t think I’ve ever seen it.
In both tests, there is a clear pattern of lower fasting BG in ABX mice. This aligns with our previous findings. We can combine the fasting BG levels at both tests (since both are baselines, just a day apart), and run an LME to test the significance of these fasting differences.
[1] “The linear mixed-effects model shows that the effect is significant. p = 1.98e-05 .” As expected, there is a significant effect of ABX-status on fasting blood glucose levels (p = 0.00002).
OGTTs Each hormone longitudinally in mice with a conventional gut
microbiome:
Are any of these differences statistically significant? Since this is longitudinal changes, I can run paired t-tests.
PGH Non ABX Baseline to Endpoint OGTT: [1] “The t-test result is not significant. p = 0.0646 .”
PL Non ABX Baseline to Endpoint OGTT: [1] “The t-test result is significant. p = 0.013 .”
Saline Non ABX Baseline to Endpoint OGTT: [1] “The t-test result is not significant. p = 0.951 .”
Are any of these differences statistically significant? Not that while for the OGTTs, I could run paired t-tests, since some ITTs failed, I can only run unpaired t-tests.
PGH Non ABX Baseline to Endpoint ITT: [1] “The t-test result is significant. p = 0.0437 .”
PL Non ABX Baseline to Endpoint ITT: [1] “The t-test result is not significant. p = 0.995 .”
Saline Non ABX Baseline to Endpoint ITT: [1] “The t-test result is not significant. p = 0.111 .”
PGH ABX Baseline to Endpoint OGTT: [1] “The t-test result is not significant. p = 0.133 .”
PL ABX Baseline to Endpoint OGTT: [1] “The t-test result is significant. p = 0.00985 .”
Saline ABX Baseline to Endpoint OGTT: [1] “The t-test result is not significant. p = 0.297 .”
PGH ABX Baseline to Endpoint ITT: [1] “The t-test result is not significant. p = 0.752 .”
PL ABX Baseline to Endpoint ITT: [1] “The t-test result is significant. p = 0.0133 .”
Saline ABX Baseline to Endpoint ITT: [1] “The t-test result is not significant. p = 0.366 .”
In each of the six groups, I want to see if there is a longitudinal change in fasting BG.
Similar to what I did for baseline, I can combine the OGTT and ITT fasting BG values and run LMEs with mouse ID as a random effect (i.e., 4 values / mouse, 2 at each time point) to increase my power to detect longitudinal changes.
Note that particularly in the ABX mice, I predict that fasting BG may decrease in all groups, including saline, because even though they are on ABX at both timepoints, there may be an additive effect of longer-term ABX exposure on BG levels.
I’m actually not sure by looking at the figures that any will be significant, but we’ll see!
PGH Non-ABX Baseline to Endpoint: [1] “The linear mixed-effects model shows that the effect is not significant. p = 0.67 .”
PL Non-ABX Baseline to Endpoint: [1] “The linear mixed-effects model shows that the effect is not significant. p = 0.304 .”
Saline Non-ABX Baseline to Endpoint: [1] “The linear mixed-effects model shows that the effect is not significant. p = 0.254 .”
PGH ABX Baseline to Endpoint: [1] “The linear mixed-effects model shows that the effect is not significant. p = 0.623 .”
PL ABX Baseline to Endpoint: [1] “The linear mixed-effects model shows that the effect is not significant. p = 0.282 .”
Saline ABX Baseline to Endpoint: [1] “The linear mixed-effects model shows that the effect is not significant. p = 0.649 .”
So NO groups saw significant changes in fasting BG from baseline to endpoint. Despite some clear changes in tolerance test responses. That’s interesting!
Endpoint OGTTs in each hormone group:
We can include all 3 hormone groups together, or to make seeing
patterns visually easier, separate each hormone to just be vs
saline
Endpoint ITTs in each hormone group:
Now we can test the stats for each of these AUC pairings.
Non-ABX
PGH vs Saline
OGTT [1] “The t-test result is not significant. p = 0.692 .”
ITT [1] “The t-test result is not significant. p = 0.627 .”
PL vs Saline
OGTT [1] “The t-test result is not significant. p = 0.465 .”
ITT [1] “The t-test result is not significant. p = 0.472 .”
ABX
PGH vs Saline
OGTT [1] “The t-test result is not significant. p = 0.803 .”
ITT [1] “The t-test result is not significant. p = 0.869 .”
PL vs Saline
OGTT [1] “The t-test result is not significant. p = 0.135 .”
ITT [1] “The t-test result is significant. p = 0.0475 .”
***
### Percent Change Responses
Now let’s see if there are statistically significant differences between the hormone groups in endpoint fasting BG values.
Non-ABX
PGH vs Saline [1] “The linear mixed-effects model shows that the effect is significant. p = 0.00257 .”
PL vs Saline [1] “The linear mixed-effects model shows that the effect is significant. p = 0.0497 .”
ABX
PGH vs Saline [1] “The linear mixed-effects model shows that the effect is significant. p = 0.0294 .”
PL vs Saline [1] “The linear mixed-effects model shows that the effect is not significant. p = 0.772 .”
Now we will test the significance for each antibiotic group pairing. I’m guessing these may all be significant.
PGH ABX vs Non ABX
OGTT [1] “The t-test result is significant. p = 0.0246 .”
ITT [1] “The t-test result is significant. p = 0.00461 .”
PL ABX vs Non ABX
OGTT [1] “The t-test result is significant. p = 0.00881 .”
ITT [1] “The t-test result is significant. p = 0.00882 .”
Saline ABX vs Non ABX
OGTT [1] “The t-test result is significant. p = 0.000758 .”
ITT [1] “The t-test result is significant. p = 0.0351 .”
As expected, for all 3 hormone groups, there is a signif dif in AUC for both OGTT and ITT between mice with and without a normal gut microbiome.
Basically, it looks like the metabolic differences between ABX and
Non-ABX groups are largely driven by changes in fasting BG, so
converting to % change makes them go away. But the metabolic differences
between hormone groups are NOT driven by fasting BG, so they largely
remain or sometimes magnify upon converting to % change.
Now let’s see if there are statistically significant differences between the abx groups for each hormone in endpoint fasting BG values. Based on what we’ve seen so far, we can expect them all to be significant.
PGH ABX vs Non ABX [1] “The linear mixed-effects model shows that the effect is significant. p = 0.000278 .”
PL ABX vs Non ABX [1] “The linear mixed-effects model shows that the effect is significant. p = 0.00312 .”
Saline ABX vs Non ABX [1] “The linear mixed-effects model shows that the effect is significant. p = 0.00287 .”
The EchoMRI produces the following measurements:
These are recorded as absolute masses in grams, but I have also converted them to percentage of body mass for an additional way to assess body composition.
Similar to the metabolic tests, each mouse was measured at baseline and endpoint. Therefore, we can assess the effects of hormone treatment via:
Longitudinal changes can be visually represented in a few dif ways:
NOTES:
Because the mice are getting older, we might expect absolute lean mass to increase across the board. That’s why looking at % can be useful too.
There are some outliers I could consider taking out. I didn’t for now…
Does the relationship between these variables and total body weight
differ by ABX status?
No clear differences in relationship by ABX status.
Let’s plot absolute lean mass (g) vs absolute body weight (g) in each group.
To see if there are differences in this relationship between hormone
groups, we will have to filter to only include the endpoint values. This
means we will be working with very low sample sizes for a regression, so
these are probably going to be not informative graphs. But nonetheless
here they are!
How about lean mass percentage vs absolute body weight (g)?
Let’s plot absolute fat mass (g) vs absolute body weight (g)
How about free mass percentage vs absolute body weight (g)?
Not that interesting.
This section is a simple look at both rates of weight gain and endpoint body mass across the six groups.
NOTE: The longitudinal visualizations below reflect actual
recorded weight. No corrections were made for cecal mass.
We can use LMEs to see if the rate of change in body mass is significantly different between different groups.
First let’s see within the Non-ABX mice if hormone affects the rate:
Though they don’t look visually distinct to me, some points are significant.
At Days 0 and 4, there is a significant difference between PGH and Saline. Day 2 is also approaching significance (p = 0.0567).
At Days 3, 4, and 5, Weight gain rates in PL and Saline are significantly different.
Next let’s see within the ABX mice if hormone affects the rate:
There is a significant difference between PGH and Saline at Days 6 and 7. Day 5 is also approaching significance (p = 0.066).
There is a significant difference between PL and Saline at Days 2, 7, and 8
Now we can see within each hormone group if the Non_ABX and ABX mice significantly differ in weight gain rates (remembering that these LMEs are not considering cecal mass).
As a simple non-interaction fixed effect, antibiotics affects weight gain in pgh mice (p = 0.001). As an interaction effect, antibiotics are significant at every day.
As a simple non-interaction fixed effect, antibiotics do not significantly affect weight gain in pl mice.
As an interaction effect, the weight gain slopes significantly differ at days 3, 4, 6, 7, and 8.
As a simple non-interaction fixed effect, antibiotics do not significantly affect weight gain in saline mice.
As an interaction effect, the weight gain slopes do not significantly differ at any days.
By D8, we can assume that cecum mass is contributing significantly to body weight in the antibiotic mice. For these endpoint comparisons, we will calculate them both corrected for cecal mass and uncorrected. I’m subtracting by full cecal mass rather than cecal contents, because the effect of abx on the mass of the cecal epithelium also seems worth correcting for to me. But I could easily switch this to cecal contents!
Because our data collection plan included regular fasts, the easiest and most useful way for us to quantify and visualize and quantify food consumption is via total food consumption, rather than daily food consumption (over 24 hour periods).
To start, I’m going to plot the total food consumption of each mouse individually with a trend line. This will help us see the extent of variation in our data and point out any possible errors.
Mouse 704 is weird and needs to be filtered out!
For the next visualization, I’m going to create a faint line for each individual mouse, plus a thicker group average.
First I’ll group just by hormone, NOT considering ABX
status:
PGH mice potentially consuming less food?
Now I’ll split the above into ABX and Non-ABX and see if that alters
paterns:
Interesting! Potentially an interaction effect between ABX and hormone
status.
Next I’ll group just by ABX, NOT considering hormone status:
No obvious difference between groups, and lots of noise.
Now let’s compare ABX and non-ABX mice within the same hormone group:
Hm! Group seeing the largest effect of ABX is saline!
NOTE: I believe in style 1, the shading around the line is showing the confidence interval, not standard error. The second format shows actual standard error.
The outputs for these models are super long, so I’m just going to summmarize. Explore the child docs for the full outputs!
Now I’m gonna test if there’s a significant difference between the slopes of these food consumption lines for each group.
Recap of results: From Day 4 onward (to Day 8), the slopes of the food consumption lines for PGH and PL both significantly differ from the slope of the Saline line, indicating a reduced rate of food consumption in the hormone-injected non-abx mice.
Recap of results: From day 5 onward (to Day 8), the rate of food consumption (slope) for PL ABX mice is significantly higher than for saline ABX mice. The rate of food consumption (slope) doesn’t significantly differ between saline and PGH ABX mice.
Recap of results: The rates of food consumption do not significantly differ in the PGH ABX and Non-ABX mice.
Recap of results: The rates of food consumption do not significantly differ in the PL ABX and Non-ABX mice. However, certain days (D3, D6, and D7) are approaching significance (0.05 < p < 0.1).
Recap of results: From Day 5 onward, the rates of food consumption are significantly higher in the Saline Non-ABX mice than in the ABX mice.
The fact that there is an ABX-based effect in the saline mice but not the PGH or PL mice is interesting!
NOTE: I’m not 100% confident that these figures are showing anything worth looking at?
I next investigated the relationship between total weight gain and total food consumption. In other words, what is the relationship between the rate at which each mouse is gaining weight and the rate at which it’s total food consumption is increasing?
First I will see if the relationship differs across hormone groups:
Next I will see if the relationship differs across antibiotic groups:
Takeaway: It looks like for the same amount of food consumed,
ABX PGH mice on average gained more weight than Non-ABX PGH mice. This
relationship doesn’t seem to hold for the PL or Saline mice, suggesting
that it’s not simply driven by cecal enlargement.
NOTE: These weights are not accounting for / subtracting cecum mass.
Verifying visual observation that the hormone group slopes don’t differ:
None of the hormone slopes significantly differ for Non-ABX with saline as the reference
For the ABX mice, PGH mice had a significantly different relationship between food consumption and weight gain than saline mice, with PGH mice gaining more weight at a higher rate (i.e., gaining more weight with the same of food consumption).
There was no significant difference between the slopes for PL and Saline ABX mice.
Does the relationship between food consumption and weight gain significantly differ between Non-ABX and ABX mice in each hormone group?
PGH:
Significantly different interaction term of antibiotics for PGH mice.
NOTE: I’m a bit confused on which model is better, since the slopes look similar, just with different intercepts, so I included both above.
PL:
Significantly different interaction term of antibiotics for PL mice.
NOTE: I’m a bit confused on which model is better, since the slopes look similar, just with different intercepts, so I included both above.
Saline:
No significant interaction term of antibiotics for saline with either model type.
I’m not showing any of the plotting here because they all looked uninteresting and the sample sizes were too low to get any good insights, but for reference, I also looked to see if total food consumption and total weight gain had significant relationships (either overall or in specific groups) with the following:
The main purpose of qPCR in MFC2 was to:
broadly verify that our antibiotic cocktail successfully suppressed gut microbial density
quantify the extent of that suppression. By what order of magnitude did ABX and non-ABX groups differ in bacterial density?
Our qPCR DF includes measurements of genomes per gram of feces. These values can differ by several orders of magnitude between the ABX and non-ABX mice, making it difficult to visualize them on the same plot. To help us better compare their differences, let’s add a new column with the log-transformed values.
Seems like approximately a 50,000 fold decrease in bacterial density
from D-5 to D5.
First I’m going to include all ABX days (D-3 to D8):
Next, to look at endpoint, we can include just days 5 and day 8, when
we expect bacterial density to be maximally suppressed in the ABX mice:
There is lower bacterial density in the PL Non-ABX mice than in the
Saline Non-ABX mice when all antibiotic hormone days are included (D3,
D5, and D8). However, this difference doesn’t reach significance when
only days 5 and 8 are included. Additionally, as we can see from the
longitudinal analyses, it appears that prior to any hormone exposure,
the mice in the PL group had naturally lower fecal bacterial loads.
Motivation: Late pregnancy is often associated with a reduction in maternal gut microbial alpha diversity, a pattern similarly observed in individuals with metabolic syndrome. This overlap has led researchers to hypothesize that gestational declines in alpha diversity may contribute to the shared metabolic alterations seen in both conditions. To investigate whether gestational changes in gut microbial alpha diversity are driven by placental hormones, we conducted a longitudinal analysis in non-pregnant mice treated with placental growth hormone (PGH) or placental lactogen (PL).
For all gut microbiome analyses, I will only be focusing on the non-abx mice. The purpose of including abx mice in the MFC2 study design was to quantify the metabolic effects of each hormone with and without a typical gut microbiome. While this design also allows us to examine the effects of the ampicillin + enrofloxacin cocktail on gut microbiome composition, that’s not the focus of this project and will not be included here.
I’ll be analyzing alpha diversity through both the Chao1 and Shannon indices. These indices are calculated differently and provide complementary insights into alpha diversity. Chao1 is a richness-based index that estimates the total number of species / OTUs in a sample, accounting for rare, lower abundance taxa. Shannon index incorporates both richness and evenness, emphasizing more abundant taxa and de-emphasizing rare taxa.
NOTE: At this stage, I opted not to include other indices, such as Simpson or Observed OTUs, that measure similar aspects of alpha diversity already captured by Chao1 and Shannon. I also didn’t use analyses less common in gut microbiome literature, like Fisher
Because we are getting signals of contamination on day 8, we can opt to
make day 5 our endpoint instead.
To me, it appears that day is the biggest driver of changes in alpha diversity, with similar patterns across time in all three groups (saline, PGH, PL) for both indices.
Excluding D8, let’s look at Baseline (pre-intervention – D-3 & D0) and Intervention (D3 & D5) to see if there are any changes to alpha diversity across the groups and whether these differences vary by group.
We can also look at just D0 to D5, rather than combining D-3 and D0
as baseline and D3 and D5 as intervention:
We also can do the same as above but comparing hormone groups at each
time point or day, rather than looking longitudinally within the same
group:
We can also look at just D0 to D5, rather than combining D-3 and D0
as baseline and D3 and D5 as intervention:
We also can do the same as above but comparing hormone groups at each
time point or day, rather than looking longitudinally within the same
group:
Takeaway: There are no consistent changes in alpha diversity driven by placental hormone exposure.
I’ll be analyzing beta diversity using four commonly used distance matrices in gut microbiome studies:
Jaccard – considers only the presence or absence of taxa (focus on rare taxa)
Bray-Curtis – incorporates relative abundances (focus on common taxa)
Unweighted UniFrac – accounts for phylogenetic distances between taxa but focuses on presence/absence
Weighted UniFrac – incorporates both phylogenetic distances and relative abundances
I’ll use these distance matrices to generate ordination plots using Principal Coordinates Analysis (PCoA) to visualize patterns of community variation over time and across hormone groups.
Since the D8s cluster together and we have reason to suspect
contamination of those samples, now I’ll filter D8 out and re-do:
Longitudinal changes in saline:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_longitud_saline_distance_matrix ~ Day, data
= jaccard_longitud_saline_data) Df SumOfSqs R2 F Pr(>F)
Model 4 1.9333 0.22765 1.8422 0.001 *** Residual 25 6.5591 0.77235
Total 29 8.4925 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in saline excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_longitud_saline_nod8_distance_matrix ~ Day,
data = jaccard_longitud_saline_nod8_data) Df SumOfSqs R2 F
Pr(>F)
Model 3 1.1531 0.18296 1.4928 0.014 * Residual 20 5.1494 0.81704
Total 23 6.3024 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PGH:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_longitud_pgh_distance_matrix ~ Day, data =
jaccard_longitud_pgh_data) Df SumOfSqs R2 F Pr(>F)
Model 4 2.3273 0.26385 2.2401 0.001 *** Residual 25 6.4934 0.73615
Total 29 8.8207 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PGH excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_longitud_pgh_nod8_distance_matrix ~ Day,
data = jaccard_longitud_pgh_nod8_data) Df SumOfSqs R2 F Pr(>F)
Model 3 1.3229 0.1996 1.6625 0.001 *** Residual 20 5.3048 0.8004
Total 23 6.6278 1.0000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PL:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_longitud_pl_distance_matrix ~ Day, data =
jaccard_longitud_pl_data) Df SumOfSqs R2 F Pr(>F)
Model 4 1.8233 0.19936 1.5563 0.002 ** Residual 25 7.3226 0.80064
Total 29 9.1459 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PL excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_longitud_pl_nod8_distance_matrix ~ Day,
data = jaccard_longitud_pl_nod8_data) Df SumOfSqs R2 F Pr(>F) Model 3
1.0817 0.1535 1.2089 0.12 Residual 20 5.9648 0.8465
Total 23 7.0464 1.0000
D3 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_phylo_d3_distance_matrix ~ Hormone, data =
jaccard_phylo_d3_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.5652 0.12177
1.0399 0.387 Residual 15 4.0765 0.87823
Total 17 4.6417 1.00000
Permutation test for adonis under reduced model Permutation: free Number
of permutations: 999
adonis2(formula = jaccard_phylo_d3_pghandsaline_distance_matrix ~
Hormone, data = jaccard_phylo_d3_pghandsaline_data) Df SumOfSqs R2 F
Pr(>F) Model 1 0.3495 0.12185 1.3875 0.111 Residual 10 2.5188
0.87815
Total 11 2.8683 1.00000
D5 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_phylo_d5_distance_matrix ~ Hormone, data =
jaccard_phylo_d5_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.4421 0.11552
0.9796 0.46 Residual 15 3.3848 0.88448
Total 17 3.8269 1.00000
D8 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = jaccard_phylo_d8_distance_matrix ~ Hormone, data =
jaccard_phylo_d8_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.5380 0.11972
1.02 0.39 Residual 15 3.9561 0.88028
Total 17 4.4942 1.00000
Since the D8s cluster together and we have reason to suspect
contamination of those samples, now I’ll filter D8 out and re-do:
Longitudinal changes in saline:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_longitud_saline_distance_matrix ~ Day, data =
bc_longitud_saline_data) Df SumOfSqs R2 F Pr(>F)
Model 4 1.1286 0.34198 3.2481 0.001 *** Residual 25 2.1716 0.65802
Total 29 3.3002 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in saline excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_longitud_saline_nod8_distance_matrix ~ Day, data
= bc_longitud_saline_nod8_data) Df SumOfSqs R2 F Pr(>F)
Model 3 0.63051 0.27611 2.5428 0.01 ** Residual 20 1.65305 0.72389
Total 23 2.28357 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PGH:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_longitud_pgh_distance_matrix ~ Day, data =
bc_longitud_pgh_data) Df SumOfSqs R2 F Pr(>F)
Model 4 1.4298 0.3745 3.742 0.001 *** Residual 25 2.3881 0.6255
Total 29 3.8178 1.0000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PGH excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_longitud_pgh_nod8_distance_matrix ~ Day, data =
bc_longitud_pgh_nod8_data) Df SumOfSqs R2 F Pr(>F)
Model 3 0.73498 0.26694 2.4276 0.004 ** Residual 20 2.01837
0.73306
Total 23 2.75335 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PL:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_longitud_pl_distance_matrix ~ Day, data =
bc_longitud_pl_data) Df SumOfSqs R2 F Pr(>F)
Model 4 1.1638 0.28962 2.5482 0.001 *** Residual 25 2.8545 0.71038
Total 29 4.0184 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PL excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_longitud_pl_nod8_distance_matrix ~ Day, data =
bc_longitud_pl_nod8_data) Df SumOfSqs R2 F Pr(>F)
Model 3 0.65641 0.21504 1.8263 0.029 * Residual 20 2.39611 0.78496
Total 23 3.05252 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
D3 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_phylo_d3_distance_matrix ~ Hormone, data =
bc_phylo_d3_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.21906 0.13519
1.1724 0.316 Residual 15 1.40136 0.86481
Total 17 1.62042 1.00000
Permutation test for adonis under reduced model Permutation: free Number
of permutations: 999
adonis2(formula = bc_phylo_d3_pghandsaline_distance_matrix ~ Hormone,
data = bc_phylo_d3_pghandsaline_data) Df SumOfSqs R2 F Pr(>F)
Model 1 0.15969 0.1705 2.0555 0.086 . Residual 10 0.77688 0.8295
Total 11 0.93656 1.0000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
D5 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_phylo_d5_distance_matrix ~ Hormone, data =
bc_phylo_d5_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.11972 0.10938
0.9211 0.513 Residual 15 0.97478 0.89062
Total 17 1.09450 1.00000
D8 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = bc_phylo_d8_distance_matrix ~ Hormone, data =
bc_phylo_d8_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.18981 0.12353
1.0571 0.387 Residual 15 1.34668 0.87647
Total 17 1.53649 1.00000
Since the D8s cluster together and we have reason to suspect
contamination of those samples, now I’ll filter D8 out and re-do:
Longitudinal changes in saline:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_longitud_saline_distance_matrix ~ Day, data
= unifrac_longitud_saline_data) Df SumOfSqs R2 F Pr(>F)
Model 4 0.09385 0.26276 2.2276 0.001 *** Residual 25 0.26332
0.73724
Total 29 0.35717 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in saline excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_longitud_saline_nod8_distance_matrix ~ Day,
data = unifrac_longitud_saline_nod8_data) Df SumOfSqs R2 F
Pr(>F)
Model 3 0.048339 0.16754 1.3418 0.057 . Residual 20 0.240176
0.83246
Total 23 0.288515 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PGH:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_longitud_pgh_distance_matrix ~ Day, data =
unifrac_longitud_pgh_data) Df SumOfSqs R2 F Pr(>F)
Model 4 0.10451 0.28253 2.4611 0.001 *** Residual 25 0.26540
0.71747
Total 29 0.36990 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PGH excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_longitud_pgh_nod8_distance_matrix ~ Day,
data = unifrac_longitud_pgh_nod8_data) Df SumOfSqs R2 F Pr(>F)
Model 3 0.049562 0.17479 1.4121 0.055 . Residual 20 0.233980
0.82521
Total 23 0.283542 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PL:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_longitud_pl_distance_matrix ~ Day, data =
unifrac_longitud_pl_data) Df SumOfSqs R2 F Pr(>F)
Model 4 0.09516 0.23233 1.8915 0.002 ** Residual 25 0.31445
0.76767
Total 29 0.40961 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PL excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_longitud_pl_nod8_distance_matrix ~ Day,
data = unifrac_longitud_pl_nod8_data) Df SumOfSqs R2 F Pr(>F) Model 3
0.04652 0.14532 1.1335 0.267 Residual 20 0.27361 0.85468
Total 23 0.32014 1.00000
D3 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_phylo_d3_distance_matrix ~ Hormone, data =
unifrac_phylo_d3_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.03878 0.1534
1.359 0.105 Residual 15 0.21402 0.8466
Total 17 0.25280 1.0000
Permutation test for adonis under reduced model Permutation: free Number
of permutations: 999
adonis2(formula = unifrac_phylo_d3_pghandsaline_distance_matrix ~
Hormone, data = unifrac_phylo_d3_pghandsaline_data) Df SumOfSqs R2 F
Pr(>F) Model 1 0.016257 0.10637 1.1904 0.282 Residual 10 0.136573
0.89363
Total 11 0.152830 1.00000
D5 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_phylo_d5_distance_matrix ~ Hormone, data =
unifrac_phylo_d5_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.028717
0.13159 1.1365 0.241 Residual 15 0.189509 0.86841
Total 17 0.218226 1.00000
D8 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = unifrac_phylo_d8_distance_matrix ~ Hormone, data =
unifrac_phylo_d8_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.015061
0.13636 1.1841 0.193 Residual 15 0.095393 0.86364
Total 17 0.110454 1.00000
Since the D8s cluster together and we have reason to suspect
contamination of those samples, now I’ll filter D8 out and re-do:
Longitudinal changes in saline:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_longitud_saline_distance_matrix ~ Day,
data = wunifrac_longitud_saline_data) Df SumOfSqs R2 F Pr(>F)
Model 4 0.14232 0.26768 2.2846 0.041 * Residual 25 0.38935 0.73232
Total 29 0.53167 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in saline excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_longitud_saline_nod8_distance_matrix ~
Day, data = wunifrac_longitud_saline_nod8_data) Df SumOfSqs R2 F
Pr(>F)
Model 3 0.12355 0.30314 2.9 0.015 * Residual 20 0.28403 0.69686
Total 23 0.40759 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PGH:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_longitud_pgh_distance_matrix ~ Day, data =
wunifrac_longitud_pgh_data) Df SumOfSqs R2 F Pr(>F)
Model 4 0.13844 0.23717 1.9431 0.078 . Residual 25 0.44527 0.76283
Total 29 0.58371 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PGH excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_longitud_pgh_nod8_distance_matrix ~ Day,
data = wunifrac_longitud_pgh_nod8_data) Df SumOfSqs R2 F Pr(>F)
Model 3 0.12596 0.25603 2.2942 0.053 . Residual 20 0.36601 0.74397
Total 23 0.49197 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PL:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_longitud_pl_distance_matrix ~ Day, data =
wunifrac_longitud_pl_data) Df SumOfSqs R2 F Pr(>F)
Model 4 0.14866 0.25557 2.1457 0.017 * Residual 25 0.43303 0.74443
Total 29 0.58169 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
Longitudinal changes in PL excluding D8:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_longitud_pl_nod8_distance_matrix ~ Day,
data = wunifrac_longitud_pl_nod8_data) Df SumOfSqs R2 F Pr(>F)
Model 3 0.14569 0.27661 2.5492 0.015 * Residual 20 0.38101 0.72339
Total 23 0.52670 1.00000
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
D3 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_phylo_d3_distance_matrix ~ Hormone, data =
wunifrac_phylo_d3_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.024242
0.09833 0.8179 0.555 Residual 15 0.222289 0.90167
Total 17 0.246531 1.00000
Permutation test for adonis under reduced model Permutation: free Number
of permutations: 999
adonis2(formula = wunifrac_phylo_d3_pghandsaline_distance_matrix ~
Hormone, data = wunifrac_phylo_d3_pghandsaline_data) Df SumOfSqs R2 F
Pr(>F) Model 1 0.015955 0.11451 1.2931 0.259 Residual 10 0.123383
0.88549
Total 11 0.139338 1.00000
D5 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_phylo_d5_distance_matrix ~ Hormone, data =
wunifrac_phylo_d5_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.011666
0.08038 0.6556 0.738 Residual 15 0.133464 0.91962
Total 17 0.145130 1.00000
D8 clustering:
STATISTICAL ANALYSIS: Permutation test for adonis under reduced model Permutation: free Number of permutations: 999
adonis2(formula = wunifrac_phylo_d8_distance_matrix ~ Hormone, data =
wunifrac_phylo_d8_data) Df SumOfSqs R2 F Pr(>F) Model 2 0.033458
0.12389 1.0606 0.385 Residual 15 0.236595 0.87611
Total 17 0.270053 1.00000
NOTE: For most distance matrices, I was able to get significant clustering by hormone group when including multiple intervention days and accounting for each individual mouse as a random effect. However, this was not resulting in R^2 values any different than when mouse ID was the only dependent variable included, so I believe it was driving the significance and creating the illusion of significant group differences. I’d love to talk through those stats methods more to ensure I’m doing everything correctly.
Key Takeaways: While there was a significant pattern of PERMANOVA separation by Day within each hormone group, indicating a longitudinal effect, there was no consistent pattern of separation by hormone group at any time point.
I can extract the data on principal coordinates from the “vectors” part of an ordination object. I can then plot these vectors across other variables, like time.
I’ll start with plotting these in terms of Bray-Curtis ordination:
Time can also be plotted against axes derived from our other distance matrices.
phyloseq Relative Abundance AnalysesIn this section, I will quantify and visualize group differences and longitudinal changes in the relative abundance of targeted taxa.
I will start with a high level look at general taxonomic groups, and then zero in on taxa in the literature that have previously been identified as differentially abundant either in murine pregnancy or in murine models of metabolic syndrome or GDM.
I will then do some visualizations of overall taxonomic composition
Changes / differences in individual taxa can be presented a number of ways. Below shows some experimentation with different visualization styles — I’m undecided on which tell the “story of the data” best.
Additionally, due to my caution around my Day 8 samples, I will sometimes opt to visualize with D5 set as endpoint.
First let’s investigate broad strokes changes in the relative abundances of two major bacterial phyla: Bacteroidetes and Firmicutes.
For our first plot, I’ll visualize relative abundance in each group
across all days
We can also show each hormone in it’s own panel
Next I’ll visualize this as baseline to intervention changes
I’m curious to visualize as a bar graph showing changes from baseline (D-3) to endpoint (D8).
To be comprehensive, I will then do the same with baseline (D-3) to
endpoint (D5), since the D8s are suspicious.
Stats testing is hidden from this html output, but these are the key
results:
No paired t-tests are significant for differences in Bacteroidetes relative abundance at D-3 and D8 for any of the hormone groups
For differences in Bacteroidetes relative abundance at D-3 and D5, the paired t-test is significant for PL mice (p = 0.036) but not for saline or PGH.
Next I’m going to calculate and plot the absolute change in Bacteroidetes % for each mouse
Change in Bacteroidetes from D-3 to D8 (i.e., D8 rel abundance -
D-3 rel abundance)
The stats analyses are hidden in this HTML, but with this group comparision in Bacteroidetes change, error bars are high enough that no two groups reach significance in difference.
Next I’m going to repeat everything that I did above, but for Firmicutes.
Stats testing is hidden from this html output, but these are the key
results:
No paired t-tests are significant for differences in Firmicutes relative abundance at D-3 and D8 for any of the hormone groups
For differences in Firmicutes relative abundance at D-3 and D5, the paired t-test is significant for PL mice (p = 0.035) but not for saline or PGH.
This matches the results found for Bacteroidetes.
Though its a bit of an outdated metric, we can also quantify the B:F
ratio:
Those have some massive error bars on the baseline days!
Stats testing is hidden from this html output, but these are the key results:
No paired t-tests are significant for differences in Bacteroidetes relative abundance at D-3 and D8 for any of the hormone groups
No paired t-tests are significant for differences in Bacteroidetes relative abundance at D-3 and D5 for any of the hormone groups, but for PL mice, it is approaching significance (p = 0.098)
To increase my sample size, I’m going to combine D-3 and D0 into
baseline and have D3 and beyond as intervention:
Stats testing is hidden from this html output, but interestingly, it seems that because the baseline error bars are so high, there is no significant difference between baseline and intervention for F:B ratio for any of the hormone groups. Potentially also because these had to be done as unpaired rather than as paired t-tests, reducing our power.
I decided to look at the following taxonomic groups based on what I’d
seen in the literature:
I’m not going to show every possible plot for each of the above 7 taxa, but I’lll show some highlights for each below: (complete rel and abs abundance analyses are found in Grace’s HEB 114 Lab 2 Draft RMarkdown
This might look like a bloom in PL mice, but when we scrutinize more
closely, it’s likely driven by a single individual (one in PL and one in
saline), rather than being a consistent trend:
From this plot, I decided that when abundances are predominantly low,
it’s better to plot each mouse individually to see if blooms are
constrained to one individual.
Seems to bloom in all groups on day 8 (day with likely
contamination):
Given the distrust of D8 samples, I’m going to filter them out:
(Note: There was one clear crazy outlier mouse, so I filtered it
out)
This also increased across groups on the day of likely contamination
But if you look at relative and absolute abundace without D8, it
might be more interesting:
It looks like it might decrease in PGH? But unclear if this is
significant at all…
This is another that spikes in absolute abundance on D8 in all
groups:
Below I’ll show relative abundance both with and without the outlier
mouse (just by changing y-axis scale)
So this taxa increases in several of the PL individuals (3), although
one mouse started with an abnormally high concentration?
Seems to have decreased in PL mice, but started from baseline at an
abnormally high level.
I’m going to start by looking at phylum-level taxonomic composition across groups.
Because we saw potential shifts in the F:B ratio of PL mice, let’s start at the phylum level:
I’ll now do the same for the levels of class, order, family, and
genus:
Maaslin2 AnalysesThis section compiles and summarizes all Maaslin2 analyses conducted for MFC2.
Rather than rerunning the analyses and regenerating output files, this document presents a comprehensive overview of the results using figures and tables generated in prior files (both on github for easy reference). This avoids redundant file creation and allows me to embed Maaslin2 outputs directly into this file.
Within each category, Maaslin2 analyses were conducted at all taxonomic levels.
Analysis categories include:
Cross-Sectional Analyses – Identifying significant microbial differences between PL or PGH groups and saline controls at a given time point
Longitudinal Analyses – Identifying significant microbial differences over time within each group
Pre-Post Intervention Analyses – Identifying microbial changes in response to hormonal intervention (comparing pre- vs. post-intervention time points, incorporating multiple measurements per mouse)
Microbiome-Metabolic Correlation Analyses – Assessing relationships between specific microbial taxa and various metabolic metrics:
Glucose metabolism (endpoint OGTT and ITT AUCs, fasting blood glucose levels)
Body composition via EchoMRI (lean, fat, total, and free water)
Weight (absolute and relative to baseline)
The above analyses were conducted both in terms of within-group variation (i.e., does individual microbiome variation correlates with individual metabolic variation within a group?) and cross-group variation (i.e., testing whether study-wide metabolic variation influences microbiome composition more than hormonal group affiliation)
Each section will include a brief methodological explanation before presenting the corresponding results.
When comparing PGH or PL mice to Saline mice including all intervention day samples (D0 and beyond), we get no significant results.
I did not test for significant microbial differences between baseline and intervention days for mice given just saline, because I wasn’t sure what the value would be at this stage. But I am happy to if we think it would be beneficial!
I tested for significant longitudinal microbial differences in mice given PGH.
Here is the statistical info on the significant result above:
Increases were seen in the phyla Tenericutes and Bacteroidetes
Decreases were seen in the phylum Actinobacteria
At the class level, nothing popped up as significant after correcting for FDR, but it appears that enrichment in Tenericutes is driven by the class Mollicutes (unadjusted p = 0.022).
At the order level, nothing popped up as significant after correcting for FDR, but it appears that enrichment in Mollicutes is driven by an order called “RF39” (unadjusted p = 0.022).
Note: Many of the family-level trends here could then be differentiated to the genus level. To make this document easier to follow, when applicable, I am going to nest those trends here, presenting the family information and then the related genus-level plots directly after.
This trend was predominantly driven by the genus Streptomyces.
Though directions of effect differed, significant effects were found in five different families within the order Clostridiales
(here’s another strong figure showing that Lachnospiraceae
decrease, but from day 0 to 5):
This trend was driven by decreases in three different genera within this family:
This trend was predominantly driven by the genus Pseudoramibacter.
This trend was predominantly driven by a genus noted as “rc4.4”.
(here’s another strong figure showing that
Dehalobacteriaceae decrease, but from day 0 to 5):
This trend was predominantly driven by the genus Dehalobacterium.
(here is also a figure for Clostridiaceae increases in PGH
mice from days 0 to 5):
This trend was predominantly driven by a genus noted as “SMB53”.
Here is the statistical info on the significant result above:
8 of the genus-level analyses were included above under their corresponding family-level analyses that also were significant.
Excluding those, there were significant decreases in 4 genera:
Here is the statistical info on the significant result above:
(Through D5): None of the groups identified above could be narrowed down by Maaslin beyond the genus level, so the species-level plots are all identical to those above.
(Through D5): We saw significant longitudinal decreases in the
following species:
We saw significant longitudinal increases in the following species:
Here is the statistical info on the significant result above:
None of these could be identified at the species level, but both decreases came from the family Lachnospiraceae, with the first plot unable to differentiate beyond family and the second able to get down to the genus Lachnospira
Increases were seen in the order Clostridiales, with Maaslin also able to distinguish that some (but not all) of these increases were driven by members of the genus Ruminococcus.
Increases were also seen in the family Lactobacillaceae.
NOTE: The weirdness of our Day 8 samples could be driving these trends. ***
I tested for significant longitudinal microbial differences in mice given PL. Because I was very skeptical of D8 values, I ran analysis at all levels that stopped at D5:
Here is the statistical info on the significant result above:
We see a mild longitudinal enrichment of Bacteroidetes and depletion of Firmicutes
We also see a mild enrichment of Verrucomicrobia
At the class level, nothing popped up as significant after correcting for FDR, but it appears that enrichment in Bacteroidetes is driven by the class Bacteroidia (unadjusted p = 0.027) and enrichment in Verrucomicrobia is driven by the class Verrucomicrobiae (unadjusted p = 0.044).
Here is the statistical info on the significant result above:
At the family level, nothing popped up as significant after correcting for FDR, but the depletion of Lactobacillales seems to be driven by the family Lactobacillaceae (unadjusted p = 0.00076).
Several other adjusted p values also were fairly low. See table below:
Significant depletion was identified in the following genera of Firmicutes:
Lactobacillus
Catenibacterium
Additionally, the following Firmicutes families showed significant depletion, but only popped up as significant in the genus-level analysis, despite not being able to be identified beyond the family level:
Lactobacillaceae (as noted above)
Erysipelotrichaceae
Here is the statistical info on the significant result above (out of order from presentation of figures above):
None of the groups identified above could be narrowed down by Maaslin beyond the genus level, so the species-level plots are all identical to those above.
To show my suspicious of the D8 samples, this is an example of what one of the species-level trends look like including those days.
To me, this trend look clearly driven by D8 contamination.
I did not test for significant microbial differences between baseline and intervention days for mice given just saline, because I wasn’t sure what the value would be at this stage. But I am happy to if we think it would be beneficial!
I tested for significant microbial differences between baseline and intervention days for mice given PGH.
Here is the statistical info on the significant result above (it’s only one result so it’s a bit boring of a table):
PGH injected-mice had a significant enrichment in Tenericutes.
Here is the statistical info on the significant result above:
The enrichment in Tenericutes is predominantly within the Mollicutes class.
Here is the statistical info on the significant result above:
The enrichment in Mollicutes is predominantly within an order called “RF 39.”
Here is the statistical info on the significant results above:
Within that order “RF 39,” we are not able to get any additional information on the family, genus, or species that is enriched. Because of this I will not show the subsequent identical figures for the lower taxonomic levels.
There is also significant depletion of a family called Peptococcaceae
Here is the statistical info on the significant results above:
Within Peptococcaceae, depletion is predominantly within a genus called “rc4.4”
Within that genus “rc4.4,” Maaslin cannot give species-level information, so the figures and stats are redundant and won’t be repeated.
I tested for significant microbial differences between baseline and intervention days for mice given PL.
Here is the statistical info on the significant results above:
Matching the results we found in phyloseq, pl-exposed mice showed an enrichment in Bacteroidetes and a depletion in Firmicutes.
Here is the statistical info on the significant result above:
It appears the depletion in Firmicutes was spread across many lower level taxa, because no classes showed significance
The enrichment in Bacteroidetes was predominantly within the class Bacteroidia
Here is the statistical info on the significant result above:
The enrichment in Bacteroidia was predominantly within the order Bacteroidiales.
At the family level, Bacteroidiales enrichment appears to have been concentrated within two families, Porphyromonadaceae and one called “S24.7.” However, while unadjusted p values were significant (p = 0.005 and 0.009, respectively), with false discovery correction, neither reached significance.
Within Porphyromonadaceae, the family that popped up as significant when p was unadjusted but not with false discovery correction, the majority of the effect came from the genus Parabacteroides (unadjusted p = 0.009)
From what I see online, Parabacteroides is considered a common gut commensal. It has also been found to “alleviate obesity and obesity-related dysfunctions in mice” (Wang et al, 2019 paper in Cell Reports), and has specifically been linked to glycemic control.
My quick lit search shows that Parabacteroides pops up a lot in GDM research, but is inconsistently found to be positively or negatively correlated. Kuang et al. (2017) did shotgun on healthy and GDM patients in the 2nd-3rd trimester (21–29 weeks) and found both that GDM patients had significantly higher Parabacteroides levels, but also that women with Parabacteroides abundance positively correlated with glucose levels during an OGTT. Their random forest model found that: “Bacterial species providing the highest discriminatory power were primarily members of the… Parabacteroides genera …consistent with our observation that Parabacteroides is the predominant genus accounting for differences in the gut microbiome between GDM patients and controls.” Dong et al. (2020) and Su et al. (2021) also found that GDM patients had a significantly higher relative abundance of Parabacteroides, as well as that “HOMA-IR increased with the higher abundance genus of Parabacteroides”. In contrast, Cortez et al. (2018) and Ma et al. (2020), found that healthy pregnant women had higher levels of Parabacteroides than GDM women. Ma et al. additionally found that women with higher Parabacteroides had lower fasting blood glucose levels.
In CD1 mice, Priyadarshini et al. (2017 – one of my fav papers) found that Parabacteroides were highest at D0 and gradually decreased during pregnancy, continuing to decrease through postpartum day 3.
Within Parabacteroides, the genus that popped up as significant when p was unadjusted but not with false discovery correction, the majority of the effect came from the species gordonii. Unadjusted p = 0.009).
I had high hopes for this but nothing popped up as significant.
The only significant effect found was at the level of order: Clostridia has a negative correlation with ITT AUC. A smaller ITT AUC indicates lower insulin resistance. So, mice with greater Clostridia levels also had greater insulin sensitivity. This effect held true when hormone both was and was not included as a random effect. The figure and stats below are for it not being included as a random effect.
Here is the statistical info for this result:
I ran these analyses two different ways: with hormone group included as a random effect, and not included.
When hormone is not included as a random effect, four taxonomic groups pop out as significant.
Here is the statistical info for these results:
When hormone is included as a random effect, the significance for those four taxonomic groups remains fairly unchanged, with the exception of the p-value for Megamonas increasing from 0.02 to 0.002.
I was unsure the ideal way to run these analyses, so I did them a few ways. First, I included lean masses from both timepoints, with individual mouse and time point included as random effects. This is basically ignoring the effects of the hormones.
NOTE: These plots all say lean mass percentage as the x axis, but are actually showing absolute lean mass. I ran these analyses before I realized the EchoMRI returns absolute and not relative lean mass values.
I only ran this analysis as the genus level (cannot remember why), and found the following:
In previous mouse studies, including Zhang et al. (2024), “SMB53 showed a stronger positive correlation with body weight, white adipose tissues, liver weight and AA metabolites, and a simultaneously negative correlation with the anti-inflammatory cytokine IL-10.”
Here is the statistical info for these results:
I then analyzed lean mass only at endpoint, accounting for hormone group as a fixed effect. The output for this was a little funky – it didn’t make plots for me, but the full results table shows some interesting things:
Here are my takeaways:
SMB53, which was positively associated with lean mass, was (unsurprisingly) also positively associated with being a PGH mouse (p = 0.035).
3 new genera popped up as having associations with lean mass:
Prevotella – positively correlated, p = 0.019
Jeotgalicoccus – positively correlated, p = 0.018 (was previously found to decrease in PGH, which is surprising)
Coprococcus – negatively correlated, p = 0.005
Similarly, I started with looking at fat masses from both timepoints, with individual mouse and time point included as random effects. This is basically ignoring the effects of the hormones.
I found a positive correlation with the Proteobacterium Gallibacterium, but it was only at non-zero levels in 6 out of 36 samples.
When I analyzed fat mass only at endpoint, accounting for hormone group as a fixed effect, I found no significant associations.
I wasn’t super interested in these, but there were tons of significant correlations with both. I’ll include the significance tables for both at the genus level below:
Free water
And here are the figures (I’ll keep them tiny):
Total water
And here are the figures (I’ll keep them tiny):
I also looked at genus-level correlations between the four body comp variables and microbial taxa, and got this heat map:
I’m not sure how to interpret it though…
I started with looking at microbial associations with weight for ALL mice, not considering hormone or day as effects, but including mouse ID as a random effect. To me, this is asking if, ignoring the experiment itself, can variation in weight across our samples be correlated with certain microbial taxa?
I ran this at the phylum and genus level, and found many significant effects for both. I won’t include ALL of the figures here, but will show the full table of significant results, and highlight some key taxa.
Phylum Level:
These results seems to correlate with commonly understood microbial correlations with adiposity: higher Firmicutes, Proteobacteria, and Actinobacteria, and lower Bacteroidetes and Verrucomicrobia.
Here is the figure for Firmicutes:
Genus Level:
Many of these taxa also popped up on our longitudinal analyses of the hormone groups (e.g., Lactobacillaceae, Faecalibacterium, Clostridium, Megamonas, etc.)
Several are also well-studied in the context of metabolism and body mass, like Akkermansia, Blautia, Bacteroides, and Lactobacillus.
Here are the plots for Lactobacillus, Akkermansia, and Bacteroides:
Lactobacillus
Akkermansia
Bacteroides
I next wanted to run analyses looking to integrate this into our actual experimental design. We can’t really parse out cause and effect, but I wanted to know if either a) Hormone exposure-driven changes in weight produce corresponding shifts in microbiome composition, or b) Hormone exposure-driven changes in microbiome composition produce corresponding shifts in weight.
I wasn’t sure what to include as fixed and random effects, and I didn’t want to miss any interesting results, so I ran a bunch of different combinations (with mouse ID always included as a random effect):
Hormone as a fixed effect with an interaction term (Weight:Hormone), and day as a random effect
Analyses run on each hormone group independently with day as a random effect
Analyses run on each hormone group independently with day as a fixed effect, filtering out the weird day 8 data
Here are the stats for analysis 1:
Phylum Level:
Genus Level:
The second analysis narrowed down our list of significant genera (see stats table below). Interestingly, one genus, “SMB53” within Clostridiaceae, which did not appear in our initial analysis list, now emerges as both positively associated with body weight AND with the PGH group. No other taxa appeared to be significantly associated with both weight and hormone group from this analysis.
Here are the associated figures for “SMB53”:
In the third analysis, day was included as a random effect and each hormone group was assessed independently.
In the PGH-dosed mice, significant weight-taxa correlations (done at genus and family level) were limited to “SMB53”, Akkermansia, Mogibacteriaceae, Bacillales, Streptophyta, and Lactobacillaceae.
In the PL-dosed mice, significant weight-taxa correlations (done at genus level) were primarily limited to members of Lactobacillales. Here is the full list:
As an example, here is the figure showing the correlation of
Lactobacillus abundance and weight in PL mice:
This is relationship is also very clear at the order level of
Lactobacillales:
In fourth analysis, day was included as fixed effect and day 8 samples were filtered out (for each hormone group independently).
In the PGH mice, “SMB53” appeared again, showing a
significant increase over time and with increasing weight:
Other than that, the PGH correlations with time and with weight largely match those discussed above.
In PL mice, there was no overlap in the taxa that significantly changed over time and that significantly explained within-group weight variation. However, this analysis also highlighted the strong longitudinal decreases in the order Clostridiales in PL mice. I re-ran the analysis at the order level to make this figure demonstrating that:
Here’s the genus-level stats:
And here’s the order-level stats: